Vortex splitting and phase separating instabilities of coreless vortices in F = 1 spinor 

Bose-Einstein condensates 
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The low lying excitations of coreless vortex states in F = 1 spinor Bose-Einstein condensates 
(BECs) are theoretically investigated using the Gross-Pitaevskii and Bogoliubov-de Gennes equa- 
tions. The spectra of the elementary excitations are calculated for different spin-spin interaction 
parameters and ratios of the number of particles in each sublevel. There exist dynamical instabilities 
of the vortex state which are suppressed by ferromagnetic interactions, and conversely, enhanced by 
antiferromagnetic interactions. In both of the spin-spin interaction regimes, we find vortex splitting 
instabilities in analogy with scalar BECs. In addition, a phase separating instability is found in the 
antiferromagnetic regime. 

PACS numbers: 03.75.Mn, 03.75.Kk, 03.75.Lm, 67.30.he 



I. INTRODUCTION 

The realization of atomic Bose-Einstein condensates 

(becs) p, a a 3 constituted the beginning of a new 

era in atomic physics. Compared with the traditional 
solid state systems, BECs in ultracold atomic gases have 
several appealing features, such as tunable interaction 
strengths, various trap potentials, and direct observation 
of the particle density. These properties offer a unique 
venue for many different types studies such as the stabil- 
ity of the BECs in a trap potential, topological defects, 
multi-component BECs, BECs in optical lattices, and low 
dimensional Bose gases [ESQ- 

In this work, we focus on vortex states in spinor 
BECs [I, i] which have been realized experimentally in 
23 Na and Rb with the hyperfme spin states F= 1 and 
F = 2 [l(| EH, E2, EE E3, UM ■ In these experiments, the 
atoms are confined optically and the condensate exhibits 
a genuine spin degree of freedom. Due to the (2F + 1) 
different hyperfme spin sublcvcls, the spinor BECs can 
have several topologically different stationary states, in- 
cluding vortex states. The quantized vortex is defined 
as a phase singularity of the condensate wave function 
[H . The phase of the wave function winds by 2irn about 
the vortex core, where the integer n is referred to as the 
vortex quantum number. For scalar BECs, the particle 
density vanishes at the vortex core due to diverging su- 
pcrfluid velocity. Due to the internal states in the spinor 
BECs, several types of vortices and other topological de- 
fects can be found. Studies related to these topologi- 
cal defects have been carried out experimentally in Rcfs. 
[HI, E3 ■ Theoretical studies of vortices and other topo- 
logical defects in F = 1 spinor BECs were initiated by 
Ohmi and Machida [|[ and Ho @. Systematic investiga- 
tions on vortices were followed by Yip [18j who considered 
both axisymmetric and nonaxisymmetric vortices, and by 
Isoshima et al. [HI, HE HH, HH , who considered only ax- 



isymmetric vortices and their excitation spectra. Studies 
of different types of related topological defects have also 
been carried out in the literature: Leonhardt and Volovik 
[23[, studied a defect referred to as Alice which is also 
known as the half quantum vortex. Stoof (24| and Mar- 
zlin et al. [2511 studied so-called skyrmions. Mizushima 
et al. j26l 1271 |2c| and Pietila et al. l29ll studied core- 
less vortices, also known as Mermin-Ho [3fJ or Anderson- 
Toulouse vortices [3l[. Furthermore, other studies of the 
exotic properties of F — 1 spinor BECs have carried out 
in Refs. [H, HE [34[. For F > 1, also many theoretical 
studies have been reported [H, H& HE HE Hi EE 511, El • 

In this work, we consider the coreless vortex states 
in spinor BECs with hyperfme spin F = 1. In the z- 
quantized basis, the condensate order parameter is de- 
noted by 4>i where i = 1,0,-1. For the coreless vor- 
tex state, the core of the vortex is filled by one of the 
components of the order parameter fa. Thus the core- 
less vortex is fundamentally different from the vortex 
in scalar BECs. Typically, the coreless vortex state 
can be defined by a combination of winding numbers 
(wi,Wo,W-\) = (0,1,2) [26|, [29[. However, by chang- 
ing the magnetization per particle M, analogous vortex 
states to the ones in a scalar BEC can be realized in 
the limit M — —1, since in this case the state </>_i is 
fully populated. In this limit, the coreless vortex state of 
the condensate corresponds to a doubly quantized vor- 
tex in a scalar condensate which is known to be dynam- 
ically unstable glllli, E3, M, US HI ■ The dy- 
namical instability is characterized by the appearance of 
the complex-frequency cigenmodes in the excitation spec- 
trum (see Sec. [TTJ). The existence of excitations with neg- 
ative but real energy is referred to as energetic instability 
or local instability and it implies that there is a stationary 
state with smaller energy to which the system tends to 
decay in the presence of dissipation. On the other hand, 
the M = 1 limit is a vortex- free state of a scalar BEC, 
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which is the ground state in nonrotating systems. It is 
thus expected that the nature of the instability of the 
corcless vortex state changes as a function of M between 
these two extreme limits. 

Let us discuss the differences between the present and 
previous studies. The condensate phase diagram in a 
plane of M and external rotation f2 has been partially 
studied in Refs. [H, H, 11]. Mizushima et al. (28|| 
focused on the ground state properties in the range 
0<M<1, and Isoshima et al. [22| studied axisymmetric 
vortex states with winding numbers Wi < 2 in the range 
— 1 < M < 1. However, these studies are not focused on 
the dynamical instability. The dynamical instability of 
the coreless vortex in a Ioffe-Pritchard magnetic field has 
been studied by Pietila et al. [29[, but only the ferro- 
magnetic case was considered. 

In this paper, we focus on the existence and character- 
istics of the dynamical instabilities in multicomponent 
systems. The coreless vortex state is an advantageous 
choice for these studies since each limit of the magneti- 
zation corresponds either to a dynamically unstable or 
stable state of a scalar BEC. We clarify how the dy- 
namical instabilities of the corcless vortex state change 
as a function of magnetization M in both ferromagnetic 
and antiferromagnetic interaction regimes. It is topical 
to study the antiferromagnetic regime since the coreless 
vortex state has been realized in 23 Na atoms with F = 1 
using the topological phase imprinting method ac- 
cording to the theoretical proposal [20l.l2l1] . On the other 
hand, 87 Rb atoms in F = 1 hyperfine spin state con- 
stitutes a ferromagnetic BEC. We demonstrate different 
aspects of dynamical instabilities in these two interaction 
regimes. The dynamical instability is suppressed by the 
ferromagnetic interactions, whereas it is enhanced by the 
antiferromagnetic interactions. In the latter case, there 
are two kinds of dynamical instabilities: the vortex split- 
ting and phase separating instabilities. In addition, we 
discuss the physical mechanisms behind these results. 

This paper is organized as follows. In Sec. [Ill we intro- 
duce a theoretical description and details of the studied 
system. In Sec. IIII1 we illustrate the condensate order 
parameter as a function of M in different spin-spin inter- 
action regimes and show a typical excitation spectrum in- 
cluding complex eigenvalues. Then we present our main 
results on the dynamical instabilities arising for different 
spin-spin interactions. In Sec. lIVi we conclude our study. 
In the Appendix, we provide a proof of the existence of 
the so-called Kohn modes in the mean-field picture of 
spinor BECs. 



field, 



H = 



dr 



i.j,k,l a 



(i) 



where 



2m 



V 2 + Vtrap(r) - n ■ (-ihr x V) - ^,(2) 



and is the bosonic field operator in the ith spin sub- 
level and m is the mass of the atoms. Here (F a ) i j is 
the (i, j) component of the spin matrix F a (a = x, y, z) 
for hyperfine spin F = 1 system. The chemical potential 
is defined as fij = fj, + jSfj, in our calculation. The sub- 
scripts {i,j,k,l} take values of the spin sublevels 1, 0, 
and — 1. The strength of the density-density and spin- 
spin interactions are denoted by the coupling constants 
g n = 47rfi, 2 (ao -(- 2a2)/3m, and g s = AkY?(o,2 — ao)/3m, 
respectively. Here ag and a2 are the s-wave scattering 
lengths between atoms with total spin and 2, respec- 
tively. In our calculation, the axisymmetric trap poten- 
tial is Vtrap(^) = ^mLd 2 r 2 with r = \Jx 2 +y 2 and the ex- 
ternal rotation is taken along the z axis f2 = (0, 0, fi). We 
consider a uniform system along the z direction. 

Following the standard procedure [1, , we write the 
time-dependent Gross-Pitaevskii (TDGP) equation as 



ift 



d 4>j(r,t) 
dt 



H? + g n J2$j(r,t)\' 



■ '1~Y,Y, (F a ) h J*(r,t)Mr,t)Mr,t). (3) 

j,k,l a 

Here, the field operator ^ has been replaced by its ex- 
pectation value tpi(r,t) = (&i(r,t)). In our simulations, 
we find the stationary state 4>i(r) using imaginary time 
propagation. 

In an axisymmetric configuration, the wave function 
can be decomposed into the amplitude and phase factor 
as 

Mr) = 4>'i{r)li{e) - Mr) cxp[z(a 4 + Wi 6)]. (4) 

Following the arguments of Isoshima et al. [22l | , station- 
ary states obey conditions 



2ao = cti + a-i + mr, 



(5) 



II. SYSTEM AND FORMULATION 



2w = wi +' 



(6) 



We begin with the second quantized Hamiltonian for 
an F = 1 spinor BEC [3] in the absence of a magnetic 



with neZ. In Eq. ([5]), we choose n = 0, ao — ot±i = 
for the ferromagnetic interaction, and n = 1, «o = t/2, 
a±i = for the antiferromagnetic case, without loss of 
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generality. Hence, we can choose {$} to be real posi- 
tive valued functions in the following discussion. These 
choices have essentially no effect in the discussion be- 
low. Note that the coreless vortex states are defined as 
(wi , wq , W- 1 ) = (0 , 1 , 2) , which satisfies Eq. ([6]) . 

Let us consider small fluctuations in the vicinity of the 
stationary state 4>: 



^(r,t) = 4 (r) + A^(r)e l - J -< ?; (r) e - t - t j(.7) 

By linearizing Eq. ([3]) with respect to A, we obtain the 
Bogoliubov-de Gennes (BdG) equation, 



£qW q 



(8) 



where the BdG operator f is composed of 3 x 3 complex 
matrices P and Q, 



T = 



P -Q 

Q* -P* 



(9) 



The matrix elements of P and Q are given by [1| 



k,l a 



Q 



kj a 



'■J 



In Eq. |(5J), the cigenfunction is denoted by 





«q,l 




Uq,l 


liq = 


Uqfi 


. «q = 


W q,0 




_ Wq,-1 _ 




_ «q -1 _ 



(12) 



Since the BdG matrix @ is generally non-Hermitian, the 
eigenvalues Eq can be complex. 

The BdG matrix has two symmetries 



and 



T* = -fiTfi 



T f = f 3 Tf 3 , 



(13) 



(14) 



where we have introduced the first and third Pauli ma- 



trices as fi = 



70 

to Q 



and r 3 = 



ro Q 

Q -T 



where 



To = diag(l, 1,1) and is a 3 x 3 matrix of zeros. The 
first symmetry in Eq. (| 13() implies the existence of two 
symmetric eigenmodes 



(Eq,Vf c 



(15) 



These modes have opposite angular momenta under the 
axial symmetry. 



Using the second symmetry in Eq. (|14[) . we obtain 



{E%-E*) J drw q (r)f 3 w q ,(r) = 0. 



(16) 



For a real eigenvalue Eq — Eq, Eq. ([T5]) implies that for 
q / q' the two modes are orthogonal. Thus we use nor- 
malization 



drw q (r)f 3 w q /(r) 



3 q,q'- 



(17) 



for quasiparticle amplitudes corresponding to real eigen- 
values of the BdG equation. Two modes provided by 
the symmetry in Eq. (|13|) give identical contribution to 
the energy of the quasiparticles and thus only the mode 
with positive norm in Eq. (|17|) is chosen as a physically 
meaningful mode. 

In the case of the complex eigenvalue E 1 * 7^ Eq, Eq. (JTHJ) 
gives J drw q f 3 Wq = and we take the following normal- 
ization condition [53j 

J drw q (r)f 3 w q <(r) 



3 qq'- 



(18) 



We use a pair of eigenmodes (Eq, w q ) and (Eq, w a ), 
for which the eigenvalues satisfy the condition i?*=_E q . 
We can find such an eigenmode as follows. By intro- 
ducing a unitary matrix U = diagL4, A*}, where A = 
e Ml diag[l, e %a , e ], and u = u$ — u\, which renders the 
BdG matrix ^ real, Eq. (jHJ) can be written in the form 



T X q Eq~X.q, 



(19) 



where T" = WTIJ is a real matrix. The complex conju- 
gate of the equation assumes the form 



f'x! 



-^q X q> 



(20) 



For a real eigenvalue £ , *=i? q , the cigenfunction can be 
taken to be real x* =x q . Hence Eqs. (|T9|) and (|20|) are 
identical. For a complex eigenvalue Eq^ Eq, eigenfunc- 
tion x q is complex, i.e., x*^x q . The eigenfunctions in 

Eq. (US) and {20} are x q = [)t Wq for Eq and x* = J7w* 
for Eq, where we used U* = . Here (_E q , w q ) is a so- 
lution of the eigenvalue equation (J8]). For a complex Eq, 
the eigenstates in Eq. ([5]) appear as a pairs (Eq, w q ) and 
(E*, w q ), where w q = t/ 2 w*. The additional factor U 2 in 
w q changes the relative phase between spin components. 

By choosing the normalization condition Eq. (|18|) for 
complex eigenmodes, one can construct a complete set 
with the complex-frequency modes [53| . The normaliza- 
tion condition in Eq. (|18p leaves the relative amplitude 
between w q and w q as well as their phase undetermined. 
In our study, we take equal amplitudes for w q and w q , 
that is, \uq ti (r)\ = |u q ,i(r)| and |u q ,j( r )l = |u q ,i(r)|. The 
physical interpretation of the quasiparticle amplitudes 
w q and w q corres pon ding to a complex eigenvalue is still 
an open question 
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FIG. 1: (color online) A set of complex- frequency modes in 
the complex plane. Four complex- frequency modes exist to- 
gether. 



The summary of the complex-frequency modes are 
shown in Fig. [TJ where we omit quantum indices in the 
figure. Modes 1 and 2 as well as modes 3 and 4 are 
linked by the symmetry in Eq. (fT3")) . Modes 1 and 3, 
and modes 2 and 4 in Fig. [T] are used to construct the 
normalization condition in Eq. (|18[) for a complex eigen- 
value -E* Eq. For the complex- frequency eigenmodes, 
the two modes which satisfy the conservation of the en- 
ergy and angular momentum are in resonance with each 
other (See Sec. QTLEfor details). 

For an axially symmetric system, all the eigenmodes 
of Eq. ([5]) can be classified with the quantum number 
qg € Z denoting the angular momentum with respect to 
the condensate. The cigenfunction is thus of the form 



' q. ! 



(r) = Uq,i(r)exp[i(q e +Wi)d], (21) 



(r) = Vq4r) exp[i{q e - Wl )d}. (22) 



We solve the BdG equation using the decomposition of 
Eqs. (|2ip and l|22[) to obtain the spectrum of the low- lying 
excitations. 

From this point on, we use dimensionless quantities. 
The energy is normalized by the trap energy Hu>, and 
the length is normalized by d = y 1 S/mw. In our study, 
we choose the density-density coupling constant to g' n = 
g n /(hujd 3 ) = 0.113, and spin-spin coupling constant to 
g' s = g s /(hwd 3 ) = ±0.001, ±0.01. The negative values 
of g s correspond to the ferromagnetic case and the pos- 
itive ones to the antiferromagnetic case. The values of 
the coupling constants g' n and g' s in the physical sys- 
tem [H, [5(| can be varied by tuning the trap frequency. 
In addition, g n can be changed by using Feshbach res- 
onances, and hence the ratio of g n and g s is also ad- 
justable. We note that a drawback in utilizing the stan- 
dard dc Feshbach resonance is that it tends to fix the 
magnetization of the cloud because of a required strong 
magnetic field. We assume an infinitely long axisymmct- 
ric system along the z axis, which renders the numeri- 
cal problem two dimensional. Alternatively, our results 
apply to pancake-shaped condensates, for which the co- 
herent dynamics in the tight direction can be neglected. 
Here, cylindrical coordinate r = (r, 9, z) is introduced 



and the integration in the two-dimensional plane is de- 
noted as f 2 Y)dr = J rdr J sm 8d9. The total number of 

the atoms N = Y^,i J2D ^ r 1^*1 2 = 1-5 x 10 3 d _1 is fixed. 
With this set of values of g' n and N, doubly quantized 
vortex states in scalar BECs have a dynamical instabil- 
ity [i^, 0, IH, |4(| . The magnetization is obtained from 

M=/ m dr(|&| 3 -|tf-i| 3 )/^ 

For low enough rotation frequencies, vortex lattices do 

not form, and hence Eq. (H|) holds. Thus the effect of the 
external rotation can be taken into account as a chemi- 
cal potential shift such that fij = fjf+jS/jf, // = /i+ Ml, 
and Sfi' = Sfx — Ml. In experiments, the magnetization 
per particle M is an observable, and hence the chemical 
potentials /ij can be treated as Lagrange multipliers in 
the calculation. Thus, the rotation cannot change the 
Gross-Pitaevskii (GP) solution under constant M. On 
the other hand, the external rotation changes the excita- 
tion spectrum by A-E q (fi)=-E q (fi)— -E q (fi = 0) = — Mlqg. 



III. RESULTS 

A. Coreless vortex states 

We study the coreless vortex states, defined by the 
combination of the phase windings of each component 
(u>i, Wo, W-i) = (0, 1, 2), for magnetization ranging from 
— 1 to 1, and for different strengths of the spin-spin in- 
teraction. In Fig. [2j we display typical spatial profiles of 
the order parameter for g' s = —0.001 and 0.001. In Fig. O 
the particle number TV; in different hyperfine spin states 
is presented as a function of M for different values of the 
spin-spin coupling constant g' s . 
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FIG. 2: (color online) The spatial profile of the order param- 
eter for g' a = -0.001 (a)-(d) and for g' s = 0.001 (e)-(h). The 
magnetization M is (a) -0.90, (b) -0.30, (c) 0.33, (d) 0.89, 
(e) -0.90, (f) -0.30, (g) 0.34, and (h) 0.90. The solid, dashed, 
and dashed-dotted lines correspond tomF = l,0, and —1 com- 
ponents, respectively. The order parameter corresponding to 
niF — O component in the antiferromagnetic case, denoted by 
the dotted line, is purely imaginary. 

According to Isoshima et al. [22| . the spin-dependent 
term of the energy density functional can be written as 



E s {r) = |{2^ 2 (r)[^ 1 (r-)±^_ 1 (r)] 2 

+ y>?(r)-4>-M 2 }. (23) 
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Here we have assumed the phase condition 7i7_i7q 2 = 
±1 which stems from the requirement that the spin- 
dependent part of the total energy is minimized. The 
upper (lower) sign corresponds to ferromagnetic (anti- 
ferromagnetic) interaction. Equation (|23[) helps to un- 
derstand the M dependence of the order parameter for 
different values of g' s . In terms of Eq. (|2"5)) , a large magni- 
tude of the spin vector is more favorable in the ferromag- 
netic case, and oppositely, in the antifcrromagnctic case, 
the spin vector tends to vanish. By comparing panels (b) 
and (f) in Fig. [5J we observe that 4>q has larger ampli- 
tude in the ferromagnetic case and therefore enhances the 
magnitude of the spin vector. Furthermore, for a broad 
range of M, Nq is finite in the ferromagnetic case whereas 
it typically vanishes for antifcrromagnetic interactions as 
shown in Fig. [31 Moreover, the fact that the mp = — 1 
component has a different winding number to mp = 1 
component explains the asymmetry of the distributions 
in Fig. El 
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FIG. 3: (color online) The ratio of the atoms in hyperfine spin 
states and total number of the atoms Ni/N as a function 
of magnetization for g' g = —0.001 (a), g' g = —0.01 (b), g' s = 
0.001 (c), and <?^=0.01 (d). The solid, dashed, and dashed- 
dotted lines correspond to ttif = 1, 0, and —1 components, 
respectively. The total number of atoms in the 2D plane is 
fixed to jV = 1.5xl0 3 d _1 . 



B. Excitation spectra and complex- frequency 
modes 

We present a typical excitation spectrum to explain 
the mechanism behind the appearance of the dynamical 
instabilities. As observed from Eq. (?]), the fluctuation 
term grows exponentially in time when some eigenvalue 
Eq is complex. This is referred to as the dynamical in- 
stability. In such case, small perturbations about the 
stationary solution of the GP equation can render it to 
decay into another state even in the absence of dissipa- 
tion. 

A typical excitation spectrum including complex- 
frequency modes is shown in Fig. [5] In this figure, the 
horizontal axis is the angular momentum quantum num- 
ber qg of the excited state and the vertical axis is the real 



_ 2 



A 

n ° * 

< 

" a i 



□ 


i ' 
> 5 i 


1 

> 

B 


A* 




mp-1 & 
m F = □ 
m F = -1 O 
complex A 



2 ' ■ ■ ■ ' ■ ^ 

-4 -2 2 

IB 



FIG. 4: (color online) The excitation spectrum for M = — 0.9 
and g' s = —0.001. The triangle at (qg =0, Re[i5 q ] = 0) corre- 
sponds to the GP solution. The labels indicate the majority 
component in the quasiparticle amplitudes m and Vi. This 
spectrum includes the complex-frequency modes denoted by 
the triangles at (q g = ±2, Re[£ q ] = ±1.82). 



part of the excitation energy Re[£q]. The eigenstate at 
qg = and i? q = corresponds to the GP solution. The 
excitation modes are labeled by the majority component 
of the excitation, that is, for the excitation with label 
mp = i. the largest amplitude of the fluctuation is given 

W m dr[K| 2 + H 2 ]. 

The excitation spectrum shown in Fig. [5] corresponds 
to g' s = —0.001 and M = —0.9. In this case, the state 
derived from the GP equation has most of the particles 
occupying the mp = — 1 component with a winding num- 
ber w— 1=2 and a small amount of mp = 1 component 
fills the core of the vortex in the <p-\ component, see 
Fig. Eta). 

The labels also illustrate the nature of the excitation 
modes. For example, in the M = —1 limit, the mp = 1, 
0, and —1 modes correspond to longitudinal spin fluc- 
tuations, transverse spin fluctuations, and density fluc- 
tuations. However, apart from this limit, the excitation 
modes are more complicated, because of the mixing be- 
tween different spin components. 

In the spontaneous dynamical excitation of the 
complex-frequency modes, conservation of the total en- 
ergy and angular momentum must be satisfied. As de- 
picted in Fig. 21 the pair of complex-frequency modes 
with (qg = -2, Re[£ q ] = -1.82) and (qg = 2, Rc[£; q ] = 
1.82) satisfies the aforementioned constraints, and hence 
the initial state with (qg = 0, Rc[£ , q ] = 0) can sponta- 
neously decay into these two states without any dissipa- 
tion. There are also additional restrictions for the ap- 
pearance of the complex-frequency modes which will be 
discussed later. We also note that external rotation does 
not affect this condition since the excitation energies are 
shifted by —Mlqg, see Sec. [TT1 

Several complex-frequency modes are found in both 
ferromagnetic and antiferromagnetic cases. Figure [5] 
presents the imaginary part of the eigenvalues as a func- 
tion of M. We find that two types of complex-frequency 
modes can appear in the coreless vortex states: (i) a pair 
of qg ~±2 modes, and (ii) a pair of qg =±1 modes. The 
former complex-frequency mode appears in the vicinity 



6 



of M = — 1 in both ferromagnetic and antifcrromagnctic 
cases as shown in Fig. [SJ The results in the M = — 1 
limit reproduce those of the doubly quantized vortex in 
a scalar BEC 0, Q IH, |H, S3- In contrast, another 
pair of complex- frequency modes with qg = ±1 emerges 
in the antiferromagnetic interaction regime. 
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FIG. 5: (color online) The absolute values of the imagi- 
nary parts of the complex-frequency eigenvalues are shown 
as a function of the magnetization M for g' s = —0.001 (a), 
£^ = -0.01 (b), ^=0.001 (c), and ^ = 0.01 (d). The q e = ±2 
and qg = ±1 modes are indicated by dashed and solid lines, 
respectively. The insets in (c) and (d) show the detailed struc- 
ture in the vicinity of M = —1 for g' g — 0.001 and g' a = 0.01. 
There are no complex-frequency modes found for M > —0.8 
in the ferromagnetic case. 




FIG. 6: (color online) The magnetization dependence of the 
excitation energies of the qg — ±2 modes for g' B — —0.001. 
The insets show the details of the spectrum where complex 
eigenenergies appear. The eigenenergy of the qe = —2 mode 
is plotted as — Re[_E q ]. The majority component in the quasi- 
particle amplitudes m and «, of the corresponding excitation 
energies is indicated by the solid, dashed, and dashed-dotted 
lines for the tuf = 1, 0, and —1 components in qg = 2 mode, 
respectively. The dots indicate qg = —2 mode, which is dom- 
inated by the tuf — — 1 component. The complex-frequency 
modes are labeled by triangles in both qg — ±2 modes. 

Figure [H] shows the excitation energies of the qg = ±2 
modes for g' a = —0.001 as a function of M. The solid, 
dashed, and dashed-dotted lines correspond to qg = 2 
modes, for which the majority components are mjr = 1, 
0, and —1, respectively. The eigenenergy of the qg = —2 
mode is plotted as — Re[£^ q ] and denoted by dots. The 



majority component for this excitation is mp = — 1. The 
complex-frequency modes appear in the regions where 
qg = 2 and qg = —2 modes overlap, due to the energy and 
angular momentum constraints. 

Let us discuss the dependence of qg = 2 modes shown in 
Fig. O These modes are classified as quadrupolc modes, 
which give rise to the two-fold rotational symmetric de- 
formation of the condensate. In spinor BECs, due to 
the multicomponent sublevels of the order parameter, 
there are three kinds of quadrupole modes: the trans- 
verse and longitudinal spin quadrupolc modes and the 
density quadrupole mode, which correspond to the three 
lowest lines around M = 1 in Fig. [HI respectively. The 
other modes with higher energy are the higher order 
quadrupole modes. The lowest density fluctuation mode 
with qg = 2 is embed at E^ = lA5Hu around M = 1, 
which is in good agreement with _E q = \p2hw derived 
within the Thomas- Fermi approximation [48| . With in- 
creasing M, since the ground state has a finite angular 
momentum (l z ) associated with the windings (0, 1, 2), the 
energy gradually shifts as Eq(M) — E^{M = 1) oc (l z ) 
p9| and stays around _E q = 1.5fko in the whole M re- 
gion. We also note that the energy of the transverse and 
longitudinal quadrupole modes, which are shown with 
solid and dashed lines near M = —1, rapidly increase 
as M decreases because of the increase of the relative 
chemical potential difference 5fi. Since the energy of the 
lowest excitations with qg = — 2 increases with M near 
M = — 1 and the resonating qg — 2 density quadrupolc 
mode remains almost constant, the complex- frequency 
modes eventually disappear, as shown in Fig. [S] The 
complex-frequency modes appear again near M = —0.9 
since the qg = —2 excitation mode finds another mode to 
pair with such that the total energy and angular momen- 
tum conservations are satisfied. 




-0.6 -0.4 
M 

FIG. 7: (color online) The magnetization dependence of the 
lowest energy excitations with qg — —2 for different values of 
g' s . The solid, dashed, dashed-dotted, and dotted lines indi- 
cate g' a — —0.01, —0.001, 0.001, and 0.01 cases, respectively. 
The shift towards positive energy increases with g' s decreas- 
ing from positive values to negative values. The inset shows 
the absolute values of the quasiparticle amplitudes {ui,Vi} 
for g a ~— 0.001 and M =—0.95. The uo and vo are neglected 
since they are vanishingly small. 
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The complex-frequency modes appear for a clearly 
wider range of values of M in the antiferromagnetic case 
compared with the ferromagnetic case [see Fig. Ufa), (b) 
and insets of panels (c) and (d)]. To explain this ten- 
dency, we consider the lowest negative energy excitation 
with qg = —2. The excitation energy increases faster 
with increasing M in the ferromagnetic case as shown 
in Fig. [7J This tendency results from the spatial pro- 
file for the qg = —2 mode. The lowest energy excita- 
tion with qg = — 2 is mainly composed of the u_i wave 
function as shown in the inset of Fig. [7J The excitation 
wave function can be generally expanded in terms of the 
qih Bcssel function J q (r) as itq,i(r) =X^i A s J qg+Wi (k s r), 
where k s = \ s /L. Here X s is the zero point of the Bessel 
function and L is the cutoff length of the system. Since 
the Bessel function behaves as J q (r) ~ r' 9 ' near r = 0, 
and qg + ui-i = 0, we have u_i oc 7'° near r = 0. Hence, 
the lowest eigenmode at qe = — 2 is the core localized 
mode and the quasiparticle amplitude u_i(r) spatially 
overlaps with 4>i, and fills the vortex core. Due to the 
coupling term — g s \4>i | 2 w_i in the BdG matrix ([5]), the 
lowest eigenvalue increases rapidly in the ferromagnetic 
regime. In addition, the slope near M = — 1 in Fig. [7] is 
steeper in the ferromagnetic case than in the antiferro- 
magnetic case. The qe = 2 modes with positive Re[-E q ] 
in resonance with lowest qg = —2 mode are less sensitive 
to changes in M as we have discussed above. Hence, the 
complex-frequency eigenmode can appear only in narrow 
magnetization regions in the case of ferromagnetic inter- 
actions. Apart from M^—l, in the ferromagnetic case, 
the coreless vortex becomes dynamically stable. 

Let us consider the difference in the density fluctua- 
tions induced by the two complex-frequency modes. The 
perturbed density profile of each component is shown in 
Fig. [5J where the first and second rows show the den- 
sity fluctuations caused by the complex-frequency modes 
with qg = ±2 and qg = ±1, respectively. From the left 
to the right column, the density of the mp = 1,0,-1 
components are shown. The mp = component of the 
qg = ±2 complex-frequency mode is neglected, because 
its amplitude is vanishingly small. 

The qg = ±2 complex-frequency mode breaks the dou- 
bly quantized vortex in the mp = — 1 component into two 
singly quantized vortices, as shown in Fig. [5Jb). This 
mechanism of dynamical instability is equivalent to the 
dynamical instability of a doubly quantized vortex in 
scalar BECs. It has also been found in the studies of 
coreless vortices induced by external magnetic fields [2!| . 
On the other hand, the qg = ±1 complex-frequency mode, 
which appears only in the antiferromagnetic regime, has 
a fundamentally different response on the condensate. 
We categorize this kind of dynamical instability mode 
as phase separation, since this mode leads to a spatial 
separation of the could into domains of a certain compo- 
nent 4>i, 4>-i or (j)Q. Although the separation is not very 
sharp, it is clearly visible in Fig. [5] Furthermore, the 
mp = l and mp = — 1 components tend to spatially over- 
lap with each other, which is attributed to the attractive 



interaction between them due to the antiferromagnetic 
interaction, as seen in Eq. ([2"5|) . 
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FIG. 8: (color online) The density profiles |*i| 2 (left col- 
umn), |^o| 2 (center column), and |^_i| 2 (right column) per- 
turbed by excitation modes. The upper row corresponds to 
the complex-frequency qg = ±2 modes for g' s = —0.001 and 
M — —0.9, the middle row to of complex-frequency qe — ±1 
modes for g' B = 0.01 and M = 0.2, and the lower row to the 
lowest real-frequency qo — — 1 mode for g' s =0.01 and M = 0.3. 
The density profile of the tuf =0 component |*Po| 2 in upper 
row is neglected since it is vanishingly small. The complex 
modes with both positive and negative qe are equally super- 
posed since the modes appear as a result of the energy and an- 
gular momentum conservation. The field of view is lOdxlOd. 
We take A = 100 to show the essential qualitative features of 
the fluctuation. 



C. Stable modes 

In addition to the dynamical instabilities, there are 
modes with real eigenvalues even if the restrictions of the 
conservation of the total energy and angular momentum 
are satisfied. For example, (i) qg — ±2 modes for g' g 

0.01 i 1 near M-— 0.85 are shown in Fig. Efa), and (ii) 
qg = ±l modes for g^O.01 near M = 0.3 in Fig. Efc) and 
(d). 

Let us first consider the case (i). The corresponding 
modes have mp = — 1 component in majority for the 
qg = —2 mode and mp = for the qg = 2 mode. In par- 
ticular, they satisfy the condition of the conservation of 
the total energy and angular momentum. Thus they can 
in principle form an excitation with complex eigenvalue, 
and in fact, this is the case for g' s = —0.01 and certain val- 
ues of M as shown in Fig. HJb) . The difference between 
these two cases can be traced back to the stationary solu- 
tion of the GP equation. According to Fig.[3]Ja) and (b), 
4>o remains negligible for M < —0.65 and g' s = —0.001, 
but 4>q is finite in the overlapping region for g' s ~ —0.01. 
Hence the existence of a finite number atoms in the cor- 
responding qg > mode can be considered as another 
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restriction for the appearance of the dynamical instabil- 
ity. 




FIG. 9: (color online) The detailed structures of the excitation 
spectra. We show qg = ±2 modes for g' s — —0.001 (a) , and 
g' s = —0.01 (b), and qg = 1 (c), and qg = — 1 (d) modes for 
g' B = 0.01. The eigenenergies for qg — — 2 and qg — — 1 are 
plotted as — Re[_E q ]. For qg = 2 and qg = 1 modes, the majority 
components of the spectrum are indicated by solid, dashed, 
and dashed-dotted lines corresponding to m,F = l, 0, and —1, 
respectively. For the qg — —2 and qg = — 1 modes, these are 
indicated by filled squares and dots corresponding to tuf — 
0, and —1 components, respectively. The complex- frequency 
modes are labeled by triangles in all cases. 

Next, we move to the case (ii). In Fig. |Htc) and (d), 
we show excitations with purely real eigenfrequencies at 
i\7 ~ 0.3 surrounded by excitations corresponding to com- 
plex eigenvalues. To understand this behavior, we con- 
sider the related density fluctuations. The density pro- 
files in the case (ii) are shown in Fig.[8]^f)-(h). We notice 
that the fluctuation leads to the precession motion of the 
(0,1,2) coreless vortex and the phase separation appear- 
ing for M = 0.2 [Fig.[Sfc)— (e)] does not occur here. Hence 
we argue that the phase separation is a characteristic fea- 
ture of this particular type of dynamical instability. 

In addition to the above discussion on the existence of 
complex-frequency modes, we note that the existence of 
dipole modes is a general feature of the excitation spec- 
trum of a harmonically trapped many-particle system. A 
generalization [5?], HH of the Kohn's theorem [59| shows 
that these center-of-mass oscillation modes should exist 
for scalar particles with energy eigenvalue i? q = hui inde- 
pendent of the interaction strength. Thus the existence 
of Kohn modes in the theory describing the system is 
typically used to check for the validity of the approxi- 
mations made. It turns out that the dipole modes have 
exactly the energy ftw in the finite-temperature Bogoli- 
ubov approximation, in which the spatial dependence of 
thermal gas component is neglected in the GP and BdG 
equations. For so-called Popov and second-order finite- 
temperature mean-field theories, the excitation ener gy i s 
very close to, although not exactly, the trap energy [601 ]. 
In Appendix A, we present a proof that Kohn modes 
with energy i? q = fkv exist for the BdG equations we uti- 
lize independent of the magnetization, density-density, or 
spin-spin interactions. 



IV. CONCLUSIONS 

We have investigated the stability of the coreless vor- 
tex states in F = 1 spinor Bosc-Einstcin condensates. 
Namely, we have calculated the low-energy excitation 
spectra in the whole range of magnetization M by solv- 
ing the Gross-Pitaevskii and the Bogoliubov-de Genncs 
equations. The complex-frequency modes, which cause 
the dynamical instabilities, have been found in both fer- 
romagnetic and antiferromagnetic cases. 

The complex-frequency modes in the ferromagnetic 
case cause the doubly quantized vortex to decay into 
a pair of singular vortices. In addition, antiferromag- 
netic interactions were found to cause phase separation 
through dynamical instability of coreless vortices. In gen- 
eral, we found that the dynamical instabilities tend to be 
suppressed by the ferromagnetic interactions and oppo- 
sitely enhanced by the antiferromagnetic ones. We also 
note that rather slow external rotation does not have an 
effect on the dynamical instabilities for a fixed magneti- 
zation. 

In addition to the conventional energy and angular mo- 
mentum conservation, we found other restrictions for the 
appearance of the dynamical instability. One such a re- 
striction for the qe > mode is the need for a consider- 
able particle number in the component of the condensate 
order parameter to be excited. Furthermore, we found 
that only certain qg < modes can resonate with other 
modes. These correspond to the vortex splitting mode 
with qg = —2 in both interaction regimes, and the phase 
separating mode with qg = — 1 in the antiferromagnetic 
regime. Due to these constraints, a dynamically stable 
coreless vortex can exist for certain magnetizations M, 
not only in the ferromagnetic case but also in the anti- 
ferromagnetic case. Our studies can be verified exper- 
imentally in fully optically trapped spinor BECs using 
present-day techniques. 
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APPENDIX A: EXISTENCE OF KOHN MODES 

Here, we show that dipole modes exist in harmonically 
trapped spinor Bosc-Einstcin condensates described by 
the employed mean-field theory. We consider a general 
system at zero temperature. The single particle Hamil- 
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tonian is defined as, 



Mr) = -^V 2 + y trap (r), 
2m 



(Al) 



where Vt rap (r) is a general three-dimensional harmonic 
trap potential, Vtrap(f) = \ m '}2u LU 'a a ' i wnere a takes 
values x, y, and z. We introduce the following creation 
and annihilation operators, 



t = J_ (oc_ _ J _d_\ 

a ~ J2\d a Ua da ) ' 



— 1 / a 



fi -2- 



(A2) 



where d a = y/h/mu! a - The introduced operators satisfy 
the bosonic commutation relation, 

[a a ,o^/] = 8 a ,a', [a a , «a] = [a£, CL^,] = 0. (A3) 

Using this notation, the single particle Hamiltonian can 
be written in the form, 



Ho(r) = ^2 hw a (a) a a a + | 



(A4) 



We denote the order parameter for an arbitrary spin 
F BEC with the (2F+l)-dimensional vector, 

*(r) = [tf F (r),tf F _i(r),--- ,*_ F (r)] T . (A5) 
The GP equation can be written in a general form, 

[H {r)To+E(r)] *(r) + A(r)**(r) = /i*(r)(A6) 



Here the (2i 7 '-|-l)-dimensional square matrices S(r) and 
A(r) are local selfenergies. A (2F-|-l)-dimensional unit 
matrix To = diag(l, • • • ,1) is also introduced. From 
Eq. jAGF we obtain a set of two equations 



H (r)- hw a r 



4*(r) 
a Q **(r) 



-{atMr)}*(r) 
{a Q E»}**(r) 



-4{A(r)}**(r) 
a Q {A*(r)*(r)} 



(A7) 



Here we introduce tq = diag(ro,ro) and a 2 x (2F+1) 
dimensional square matrix Hq 



H (r) = diag {if (r)-/i}T + £H 



-{i?o(r)-/i}ro-£*(r) 



(A8) 



From Eq. (|A5[) one can derive the general form of the 
BdG equation 



Hn 



u v (r) 




_ v v (r) _ 


+ 



A(r)v„(r) 
-A*{r)u v (r) 



= E V 


u v (r) 




_ v v (r) _ 



.(A9) 



At zero temperature, A(r) is equal to A(r) in the GP 
equation. Let us take an ansatz 



(A10) 



and write the BdG equation using Eq. (|A7[) in the form 



u v=a (r) 




' a* ¥(r) " 


_ v v = a {r) _ 







[E a - huj a ] 







_a a <S>*{r)_ 





-{ fl tS(r)}*(r)-flt{A(r)**(r)} + A(r)« a r(r) 
{a Q S*(r)}**(r) + a Q {A*(r)*(r)} - A* (r)4 *( P ) 



(All) 



All results above are for a general BEC with hy- 
pcrfinc spin F. Below, we restrict the discussion to 
F = 1 case since the selfenergy for this case is known. 
Here, the order parameter takes the form <&(r) = 
\&i (r), ^>o(r), , 5_i(r)] T . Using the following notation 

- - 1 F 1 ' for i/= 1,2,3 ' lA12j 



9v ~ {S for £=1,2,3 ' (A13) 
the selfenergies are written as, 

E(r) = g v [^(r)A v ^(r)A v +A v ^(r)^A v ](AU) 

A(r) = -g v A v ¥(r) [&(r)A v ]* , (A15) 



where summation over repeated superscripts is implied. 
We substitute these selfenergies to the BdG equation 
([ATTj) . and using the condition ['J/t^vf] * = vpt^"^/ we 
finally observe that 



[E a - huj a ] 



' ot*(r) " 




?7(r) 









(A16) 



where we have defined, 

r,(r) = -{4E (r)Wr) _ a t { A(r)*>)} 
+ A(r)a a **(r). 

Assuming that F Q is real, Eq. (|A16j) yields 

(£ a -fc t ; a )(a],+a a )¥(r)=0. (A17) 

Since {a) a + a a )' < b{r) ^ 0, we conclude that there always 
exists a mode with energy E a — hw a . Therefore the Kohn 
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mode exists irrespective of the atom-atom interactions. numerical calculations, we typically find the dipole mode 
The eigenvector and eigenenergy are given by E a = hui a with a relative error is less than 1.5 x 10 -5 . 
and [u v ,v v ] T = [aj, 1 4'(r), a Q 1 4'*(r)] T , respectively In our 
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